//Opt_Kernel kmeans.
//Esta es una version mejorada de Kernel Kmeans
//La idea es hacer mas rapido el agrupamiento.
//Primero se divide aleatoriamente todos los puntos en los N_cent y se almacenan en partitions.
//Luego se calcula por cada cluster el centroide y  el dato mas cercano a ese centroide dentro de
//ese cluster es el nuevo centroid (Es el indice y se almacena en ind_cent).
//Finalmente, se reasignan todos los puntos en new_partitions y new_counts.
//Yo compare los resultados con Kernel Kmeans y son similares y este funciona mas rapido cuando
//se tiene a gran cantidad de datos.



inline
op_kkmeans::op_kkmeans(const mat &in_Kn, const uword in_Ncent)
: K(in_Kn), N_cent(in_Ncent), N_data(K.n_cols)
{
}

inline
void
op_kkmeans::calc(const uword max_iter)
{
  
  uword iter = 0;
  umat comparing;
  init_partitions();
  
  //only to test:
  //update_kkmeans();
  //partitions = new_partitions;
  
  //while ( err > 0.005 && iter < max_iter)
  do 
  {  
    cout << "Doing iter= "  << iter << endl;
    update_kkmeans();
    comparing = (partitions != new_partitions);
    //partitions.print("partitions");
    //comparing.print("comparing: ");
    if (iter < max_iter)
    {
      err = sum(sum(comparing));
      
    }
    else
    {
     err = 0;
     distortion = 0; 
    }
    
    cout << "err= " << err << endl;
    cout << "distortion= " << distortion << endl;
    //err = 1;
    partitions = new_partitions;
    counts = new_counts;
    iter++;
    //cout << "(err !=0 )  || (iter < max_iter) = " << ((err !=0 )  || (iter < max_iter)) << endl;
  }while ( (distortion >0 )  ); // while ( (err !=0 )  );
  
  
  
  //cout << "Saving Lables: " << endl;
  //save_labels("labels_grass.dat");
  
}



inline
void
op_kkmeans::init_partitions()
{
  cout << "Ndata: " << N_data << endl;  
  cout << "Ncent: " << N_cent << endl;
  
  vec tmp = randu<vec>(N_data ) * (N_cent);
  //tmp.print("tmp:");
  uvec tmp2 = conv_to<uvec>::from( tmp );
  //tmp2.print("tmp2:");
  
  partitions.zeros(N_data,N_cent);
  new_partitions.zeros(N_data,N_cent);
  
  counts.zeros(N_cent);
  
  for (uword c = 0 ; c < N_cent; ++c)
  {
    uvec q1 = find (tmp2 == c);
    counts(c) = q1.n_elem;
    //q1.print("q1 bf: ");
    q1.resize(N_data);
    //q1.print("q1 after: ");
    partitions.col( c ) = q1;
  }
  //partitions.print("initial_indices:");
  //counts.t().print("counts:");
}


inline
void
op_kkmeans::update_kkmeans()
{
  
  double dist_cent;
  double a;
  double b;
  
  uword best_ci;
  vec ind_cent;
  ind_cent.zeros(N_cent);
  new_counts.zeros(N_cent);
  new_partitions.zeros(N_data,N_cent);
  distortion = 0;
  //cout << "Initialising distortion= " << distortion << endl;

  //partitions.print("initial partitions:");
  //K.print("K:");
  //counts.print("counts");
  for (uword c = 0; c < N_cent; ++c) //N_cent
  {
    double sum2 = 0;
    
    //cout << "counts(c): " << counts(c) << endl;
    vec dist(counts(c));
    
    if (counts(c)> 0)
    {
      for (uword i = 0; i < counts(c); ++i)
      {
	uword xi = partitions(i,c); 
      for (uword j = 0; j < counts(c); ++j)
	{
	  uword xj = partitions(j,c); 
	  sum2+= K( xi, xj );
	  //cout << "K( a , xi ): " << K( a , xi ) << endl;
	}
      }
	
	
      for (uword i = 0; i < counts(c); ++i)
      {
	uword a = partitions(i,c); 
	double sum1 = 0;
	
	for (uword j = 0; j < counts(c); ++j)
	{
	  uword xi = partitions(j,c); 
	  sum1+= K( a , xi );
	  //cout << "K( a , xi ): " << K( a , xi ) << endl;
	}
	//cout << "sum1= " << sum1 << " & sum2= " << sum2 << endl;
	dist_cent = K(a,a) -2*sum1/counts(c) + sum2 / pow(counts(c),2 );
	dist(i) = dist_cent;   
      }
      
      //dist.print("dist:");
      double mind = dist.min(best_ci);
      distortion+=mind;
      ind_cent(c) = partitions(best_ci,c); ///Que pasa cuando no hay datos asociados a ese centroide??
      //cout << "ind_cent(" << c << ")=" << ind_cent(c) << endl; 
    }
  }
  //ind_cent.print("ind_cent");
  
  for (uword ni = 0; ni < N_data; ++ni)
  {
    vec dist(N_cent);
    a = ni;
    
    for (uword ci = 0; ci < N_cent; ++ci)
    {
      b = ind_cent(ci);
      dist(ci) = K(a,a) - 2*K(a,b) + K(b,b);
    } 
    dist.min(best_ci);
    new_partitions( new_counts(best_ci) , best_ci) = ni;
    new_counts(best_ci) += 1;
    
    /* new_partitions(new_counts(best_ci),best_ci) = ni;
    *    new_counts(best_ci) += 1;*/
  }
  
}

//This function is only useful to plot
//the labels in Matlab.
inline
void
op_kkmeans::save_labels(std::string name)
{
  uvec labels;
  labels.zeros(N_data);
  // uvec q1 = find (tmp2 == c);
  
  for (uword ci = 0; ci < N_cent; ++ci)  
  {
    uvec v = partitions.col(ci);
    //v.print("v: ");
    v.resize(counts(ci));
    //v.print("v: ");
    //labels(v).fill(ci);
    for (uword n = 0; n< counts(ci); ++n)
    {
      labels(v(n)) = ci;
    }
    
    //labels.print("labels: ");
    //getchar();*/
  }
  //labels.t().print("labels: ");
  labels.save(name, raw_ascii);
}


inline
void
op_kkmeans::print_labels(const field<std::string>  labels)
{
  //partitions.print("Final Partitions");
  field<std::string> lab_partitions( N_cent, N_data  );
  // uvec q1 = find (tmp2 == c);
  
  for (uword ci = 0; ci < N_cent; ++ci)  
  {
    uvec v = partitions.col(ci);
    //v.print("v: ");
    v.resize(counts(ci));
    //v.print("v: ");
    //labels(v).fill(ci);
    cout << "Cluster " << ci + 1 << " has  " << counts(ci) << " points" << endl;
    for (uword n = 0; n< counts(ci); ++n)
    {
      //labels(v(n),0) = lab(ci,0);
      cout <<  labels( partitions( n,ci ),0 ) << endl;
      lab_partitions( ci, n ) = labels( partitions( n,ci ) ); 
      // getchar();
    }
    cout << endl;
    
    //getchar();*/
  }
  //lab_partitions.print("Final labels: ");
  
}

inline
uvec
op_kkmeans::get_id_summ()
{
  uvec indices = sort_index(counts);
  uvec indices_summary;
  
  indices_summary.zeros(N_data);
  
  //partitions.print("Que pasa??");
  //indices.t().print("indices:");
  int sum = 0;
  int p_sum = 0;
  int ii   = 0;
  uword pos2;
  uword pos;
  double length_summ = N_data*30/100;
  while (sum < length_summ )
  {
    pos = indices(ii);
    if (counts(pos) > 0)
    {
    uvec tmp = partitions( span(0,counts(pos) - 1 ) , pos );
    //tmp.t().print("tmp");
    indices_summary.subvec(sum, sum+counts(pos) - 1 ) = tmp; 
    //indices_summary.t().print("indices_summary");
    ii++;
    p_sum+= counts(pos);
    pos2 = indices(ii); 
    sum+=counts(pos) + counts(pos2);
    }
    else
    {
      ii++;
    }
   
    //cout << "sum: "<< sum << endl;
    //indices_summary.t().print("Que pasa??");
  }
 
    //cout << "p_sum " <<   p_sum << endl;
    indices_summary.resize(p_sum);

  //indices_summary.t().print("indices_summary");
  return sort( indices_summary ); // returning only the smaller clusters 
}


inline
umat
op_kkmeans::get_partitions()
{
  //partitions.print("partitions");
  return partitions;
  
}

